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ABSTRACT 

The stability of a recently proposed general relativistic model of galaxies is studied 
in some detail. This model is a general relativistic version of the well known Miyamoto- 
Nagai model that represents well a thick galactic disk. The stability of the disk is 
investigated under a general first order perturbation keeping the spacetime metric 
frozen (no gravitational radiation is taken into account). We find that the stability 
is associated with the thickness of the disk. We have that flat galaxies have more 
not-stable modes than the thick ones i.e., flat galaxies have a tendency to form more 
complex structures like rings, bars and spiral arms. 
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1 INTRODUCTION 

The natural shape of an isolated self-gravitating fluid is 
axially symmetric. For this reason, exact axial symmet- 
ric solutions of Einstein field equations are good can- 
didates to model astrophysical bodies in General Rela- 
tivity. In the last decades, several exact solutions were 
studied as possible galactic m odels. Static thin d i sk so - 
lutio ns were first studi e d by uBonnor fc Sackfieldl (|l968l ) 
and iMorgan fe Morgan! l| 19691 ). where they considered 
disks without radial pressure. Disks with radial pres- 
sure and with radial t ensio n had been considered b y 
IMorgan fc Morgan! ||l970h and IGonzalez fc Letelierl l|l999h , 
respectively. Self-simila r stat ic dis ks wer e stu died by 
iLvnden-Bell fc Pineaultl l|l97Sl) . and iLemos! (|l989h . More- 
over, solutions that involve super positions of black holes 
with s tatic disks were a n alyzed by ILemos" fc Letelierl l|l993l . 
1 1994 Il996h and iKleinl (|l997l V Also, relativistic counter- 
rotating thin disks as sources o f the Kerr type metrics 
were found by iBicak fc Ledvinkal (|l993l ). Counter-rotating 
models with radial pressu re and dust disks w it hout radial 
pressure were stud ied _by IGonzalez fc Espitial ([20031 ) . and 
iGarcia fc Gonz alez ( 2004) , resp ectively; while rotating disk s 
with heat flow were studied by IGonzalez fc Letelierl |2000). 
Furthermore, static thin disks as sources of k nown vacuum 
spacetimes fr om the Chazy-Curzon metric (|Chazvl 1 1924 



Curzonl \l92$ ) and Zippy- Voor hees llZipovl Il96q; IVoo rhce; 



1970) metric were obtained b y Bicak, Lvndcn- Bell fc Katz 
|l993h . Also. lBicak. Lvnden-Bell fc Pichonl l| 19931) found~an 
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infinite number of new relativistic static solutions that cor- 
respond to the classical galactic disk p otentials of Kuzmin & 
Toomre (iKuzminll 19561; |Toomrelll963T ) and Mestel & Kalnajs 
|Mestellll963l ; lKalnaislll972l). Stationary disk mod els includ- 
ing ele ctric fields ([Ledvinka. Zofka fc Bi cak 1999), magnetic 
fields (|Letelierlll999 | ). and both el ectric and magnetic fields 
l|Katz. Bicak fc Lvnden-Belllll999l ) had been studied. In the 
last years, exact solutions for thin disks made with sin - 
gle and comp osite halos of matter |Vogt fc Letelierl 120031 ). 
charge d dust ( Vogt fc Letelierl l2004al ) and charged perfect 
fluid (IVogt fc Letelierl l2004bl ) were obtained. For a sur- 
vey on relativistic gravitating disk s, see ISemera 2 (|2002l ) 
and H aras, Hure fc Semerakl (|2004l ). Most of the models 
constructed above were found using the metric to cal- 
culate its energy momentum-tensor, i.e. an inverse prob- 
lem. Several exact disk solutions were found using the di- 
rect method that consists in computing the metric for 
a given energy momentum tensor representing the disk 
jNeugebauer fc Meinel 1 19951; iKlein fc Richtei] 1 19991 ; IKleinl 
l200ll ; lFrauendiener fc Kleinll200ll : lKleinll2002l . l2003al lbh ■ In a 
first approximation, the galaxies can be thought to be thin, 
what usually simplifies the analysis and provides very useful 
information. But, in order to model real physical galaxies 
the thickness of the disks must be considered. Exact axi- 
ally symmetric relativisti c thick disks in different coo rdinate 
systems were studied by Gonzale z fc Letelierl l|2004 ). Also, 
different thick disks were obtained from the Schwarzschild 
metric in different coordinate s systems with the "di splace, 
cut, fill, and reflect" method (IVogt fc Letelierll2005al ). 

The applicability of these disks models to any struc- 
ture found in Nature lays in its stability. The study of the 
stability, analytically or numerically, is vital to the accep- 
tance of a particular model. Also, the study of different 
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types of perturbations, when applied to these models, might 
give an insight on the formation of bars, rings or differ- 
ent stellar patterns. Moreover, a perturbation can cause the 
collapse of a stable object with the posterior appearance 
of a different kind of structure. An analytical treatment of 
the stability of disk s in N e wtonian theory can be found in 
iBinnev fc Tremaind (|l987h . iFridman fc Polvachenkol (|l984T l 
and references therein. In general, the stability of disks 
in General Relativity is done in two ways. One way is to 
study the stability of the particle orbits al o ng ge odesies. 
This kind of study was made by lLetelierl j2003) trans- 
form i ng the Rayleigh criterion of stability (|Lord Ravleighl 
Il9ld : lLandau fc Lifshit j 1 19871 ) into a general relativistic 
formulation. Using this criterion, the stability of orbits 
around black holes surroun ded by disks, rings and multi- 
polar fields were analyze d jLetelierl 120031 ) . Also , this cri- 
terion was employed by IVogt fc Letelierl (|2003h to study 
the stability of the isotropic Schwarzschild thin disk, and 
thin disks of single and composite halos. The stability of 
circular orb its in sta t ionary axis ymmetric spacetimes was 
studie d by iBardeenl (Il970l ) and |Ab ramowicz fc Prasannal 
(1990). Moreover, the stability of circular or bits of the 
Lemos-Letelier solution dLemos fc Letelierl H994T I for the su- 
pe rposition of a blac k hole a nd a fl at ring was con sidered 
by ISemerak fc Zacekl feOOOallbh an d ISemerakl i|2003l ). Also, 
iBicak. Lvnden-Bell fc Katzl l|l993l ) analyzed the stability of 
several thin disks without radial pressure or tension study- 
ing their velocity curves and specific angular momentum. 
Another way of studying the stability of disks is perturbing 
its energy momentum tensor. This way is more complete 
than the analysis of particle motions along geodesies, be- 
cause we are taking into account the collective behavior of 
the particles. However, there are few studies in the litera- 
ture performing this kind of perturbation. A general stability 
study of a relativistic fluid, with bo th bulk and dynamical 
viscosity, was done bv lSeguir] (|l975l ). He considered the co- 
efficients of the perturbed variables as constants, i.e. local 
perturbations. Usually, this condition is too restrictive. Sta- 
bility analysis of thin disks from the Schwarzschild metric, 
the Chazy-Curzon metric and Zipoy-Voorhees metric, per- 
turbing their energy momentum t ensor with a gene r al firs t 
order perturbation, were made bv lUievic fc Letelierl l|2004h . 
finding that the thin disks without radial pressure are not 
stable. Moreover, stability analysis of the static isotropic 
Schwarzschild thick disk as well as the ge neral perturbation 
equat ions for thick disks were studied bv lUievic fc Letelierl 
l|2007h . 

In Newtonian gravity, models for globular c l usters 
and spher i cal ga laxies were developed by IPlummerl (jl91ll ) 
and iKind (|l966h . In the case of disk galaxies, important 
thick disk models were obtained by Miyamoto and Nagai 
ijMivamoto fc Naeail Il975l : iNaeai fc Miyamoto! Il976l ) from 
the prior work of iKuzmini l| 19561) and iToomrd (I1963T ) about 
thin disks galaxies. Miyamoto and Nagai "thickened-up" 
Toomre's series of disk models and obtained pairs of three- 
dimen sional potential and density functions. Also, ISatohl 
(1980) obtained a family of three-dimensional axisymmet- 
ric mass distribution from the higher order Plummer mod- 
els. The Miyamoto-Nagai potential shares many of the im- 
portant properties of actual galaxies, especially the contour 
plots of the mass distribution which are qualitatively similar 
to the light distribution of disk galaxies l|Binnev fc Tremaind 



Il987l) . Recently, two different extensions of the Miyamoto- 
Nagai pot ential appeared in the literature: a triaxial gener- 
alization l|Vogt fc Le telier 2007) which has as a particular 
case the original axially symme tric model, and a relativistic 
version |Vogt fc Letelierl l2005bl ) which has as a Newtonian 
limit the same original model. 

In order to have a general relativistic physical model 
for galaxies, we must consider, first of all, the thickness of 
the disk and its stability under perturbations of the fluid 
quantities. The purpose of this work is to study numeri- 
cally the stability of the general relativistic Miyamoto-Nagai 
disk under a general first order perturbation. The perturba- 
tion is done in the temporal, radial, axial and azimuthal 
components of the quantities involved in the energy mo- 
mentum tensor of the f luid. In the general thick disk case 
|Uievic fc Letelier 2007]), the number of unknowns is larger 
than the number of equations. This opens the possibility of 
performing several types of combinations of the perturbed 
quantities. In this manuscript we search for perturbations 
in which a perturbation in a given direction of the pres- 
sure creates a perturbation in the same direction of the four 
velocity. The energy momentum perturbation considered in 
this manuscript is treated as "test matter", so it does not 
modified the background metric obtained from the solution 
of Einstein equations. 

The article is organized as follows. In Sec. [2] we present 
the general perturbed conservation equations for the thick 
disk case. The energy momentum tensor is considered di- 
agonal with all its elements different from zero. Also, in 
particular, we discuss the perturbations that will be con- 
sidered in some detail in the next sections of this work. In 
Sec. [3] we present the thick disk model whose stability is 
analyzed, i.e. the general relativistic Miyamoto-Nagai disk. 
The form of its energy density and pressures, as well as, the 
restrictions that the thermodynamic quantities must obey 
to satisfy the strong, weak and dominant energy conditions 
are shown. Later, in Sec. [4] we perform the perturbations to 
the general relativistic Miyamoto-Nagai disk; in particular 
we study its stability. Finally, in Sec. [S] we summarize our 
results. 



2 PERTURBED EQUATIONS 

The thick disk considered is a particular case of the general 
static-axially-symmetric metric 



2*i ,,2 . 2*2 r>2 7 fl 2 . 2<S 3/ , n 2 , , 2 

-e L dt +e 1 R dO + e ^{dR + dz 



(1) 



where $1, ^2 and ^3 are functions of the variables (R,z). 
(Our conventions are: G = c = 1, metric signature +2, par- 
tial and covariant derivatives with respect to the coordinate 
x 11 denoted by , \i and ; \i, respectively.) 

In its rest frame, the energy momentum tensor of the 
fluid T M " is diagonal with components (-p,pn,pe,p z ), where 
p is the total energy density and (pn,pe,Pz) are the radial, 
azimuthal and axial pressures or tensions, respectively. So, 
in this frame of reference, the energy momentum tensor can 
be written as 



pU»U v + p R X»X v + Pe Y»Y v + pzZ^Z" 



(2) 



where (7 M , A M , Y" M , and Z M are the four vectors of the or- 
thonormal tetrad 
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W = e-** (1,0, 0,0), 

X» = e-^(0, 1,0,0), 



yM _ 



-(0,0,1,0), 



^ = e -* 3 (o, 0,0,1), 



(3) 



which satisfy the orthonormal relations. Note that with the 
above definitions, the timelike four velocity of the fluid is 17 M 
and the quantities , Y^, and Z^ are the spacelike princi- 
pal directions of the fluid. Furthermore, the energy momen- 
tum tensor satisfies Einstein field equations, G^u = SttT^. 
Moreover, the quantities involved in the energy momentum 
tensor and the coefficients of the perturbed conservation 
equations are functions of the coordinates (R, z) only. Let 
us consider a general perturbation A p of a quantity of 
the form 



A P (t,R,6,z) = A"-{R,z) +5A t *(R,z)e 



i(ko 9 — wt) 



(i) 



where A^(R, z) is the unperturbed quantity and 
SA^{R,z)e l{ke<) - wt) is the perturbation. Replacing © 
for each quantity in the energy momentum tensor ((2)1 and 
calculating the perturbed energy momentum equations, 
6T£" = 0, we obtain 
p — t 

SU^ipU 1 + €,ip R X R ) + 5U z z (pU l + tzp z Z z ) 
+8U R [F(t, R, pU*) + ti, RPR X R + CiF(t, R,p R X R )] 
+5U {ik e {pU t +&p e Y e )] 
+5U z [F(t, z, pU*) + &, z p z Z z + &F(i, z, Pz Z z )] 
+Sp(-iwU t U t ) = 0, (5) 

p = R 

5 Pr m{X r X r ) + SU^-iwipU* + ZiprX r )] 

+s P (u t u t r R ) + 8p R G{R, r, x r x r ) 

+5p g (Y e Y e Tf g ) + 5 Pz (Z z Z z T R ) = 0, (6) 
p = 8 

SU^-wipU* + £,2PeY B )] + 5pg{kgY e Y B ) = 0, (7) 
p = z 

5p z , z {Z z Z z ) + SU^-iwipU* + izp z Z z )\ 



+Sp{U t U t Y z tt ) + 5p R (X R X R T z RR ) 
+Sp e (Y e Y T z eg ) + 8p z G(z, z, Z Z Z Z ) = 0. 



where 



F{I,J,K) = K,j + K(2T I IJ + r^j 
G(I,J,K) = K,j + K(Tjj + r^j) 



(8) 

(9) 
(10) 



and T§ 7 are the Christoffel symbols. In finding Eqs. ©- 
([8]) we assumed that the perturbed energy momentum ten- 
sor does not modify the background metric. Also, we disre- 
gard terms of order greater o r equal to S 2 . For details see 
lUievic fc Letelierl (|2004 l2007n . 

Besides the four equations furnished by the energy mo- 
mentum conservation equations, T.£" = 0, there is another 
important conservation equation, the equation of continuity, 



where n is the proper number density of particles. The 
proper number density of particles n, and the total energy 
density p are related through the relation, 



nm b + e, 



(12) 



where m/, is the constant mean baryon mass and e the inter- 
nal energy density. Multiplying Eq. (|12[) by performing 
the covariant derivative (; p) and using Eq. {TTJ, we obtain 
that 



(13) 



Now, from the relation U V T^ = and the energy momen- 
tum tensor ((2]), we obtain an expression for (pC/ M ) ;M . Sub- 
stituting this last expression into Eq. (|13[l we finally arrive 
to 

(d7 M ) ;fl = p R X"U v X^+peY"U v Y^+p z Z^U v Z^, (14) 

which is a first order differential equation for e. Therefore, 
with e given by (| 14f) the equation of continuity ((TTJ is satis- 
fied. For this reason, the continuity equation can be omitted 
in our stability analysis because, in principle, we can always 
find a solution for e. Hereafter, the contribution of nm b and e 
to the total energy density are taken into account in p. In the 
case in which the internal energy density of the fluid is given, 
the equation of continuity must be considered. The thermo- 
dynamic properties of the system can be obtained from ob- 
servations or theoretically, e.g. from the Fokker-Planck equa- 
tion, where we obtain the particle distribution function of 
the disk. Solving the three dimensional Fokker-Planck equa- 
tion is not an easy tas k, but some progress in New tonian 
gravity had been done (|Uievic fc Letelierl [20051 . 120061 ). 

The four equations, (J5j)- ((HJ) , contain seven independent 
unknowns, say SU R ,5U ,SU Z , 5p,5pR, 8pg,5p z . So, at this 
point, the number of unknowns are greater than the num- 
ber of equations. This opens the possibility to perform dif- 
ferent kind of perturbations. In this article we are inter- 
ested in perturbations in which the velocity perturbation 
in a certain direction leads to a pressure perturbation in 
the same direction. For example, if we perturbed the axial 
component of the velocity, SU Z , then we must perturb 5p z . 
With the above criterion, and without imposing any extra 
conditions to the individual perturbations, only four pertur- 
bations combinations are allowed and will be considered in 
our thick disk model. Furthermore, we perform the pertur- 
bation 5U R , Spr, 5U Z ,5p z with the extra imposed condition 
8pR = 5p z . In this particular case, the system of equations 
reduces to a second order partial differential equation. 



3 GENERAL RELATIVISTIC 

MIYAMOTO-NAGAI GALAXIES 

A static general relativi stic version of the Miya moto-Nagai 
disk was constructed bv lVogt fc Letelierl l|2005bl ) by making 
a correspondence between the general isotropic line element 
in cylindrical coordinates and the Miyamoto-Naga i model 
i|Mivamoto fc Nagai|[l975l ; iBinnev fc Tremainell 19871 ). These 
general relativistic disks are obtained with ((TJ and the spe- 
cializations, 



(nt/") ;M = 0, 



(11) 



*i = In 



1-/ 
1 + / 



(15) 
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$ 2 

where / 



*3 



21n(l + /), 



(16) 

m is the mass of the disk, 



2\/ R 2 +(a + yJz2+bl) 

and (a, 6) are constants that control the shape of the density 
curves. With this metric, the energy density and pressures 
for the general relativistic Miyamoto-Nagai disk are 



b 2 [all 2 + (a + £) 2 (a + 3£)] 



4tt£ 3 [|+x] 5 



Pr = pe 



b 2 [aR 2 + (a + 2 (a + 2£)] 
b 2 (a + 2 



(17) 
(18) 
(19) 



where £ = sj z 2 + b 2 and \ = \/R 2 + (a + £) 2 - Without los- 
ing generality we set m = 1 in Eqs. (|17|l - (|19| l. To satisfy 
the strong energy condition (gravitational attractive mat- 
ter) we must have that the "effective Newtonian density" 
A — p + pr + pe + p z ^ 0. The weak energy condition re- 
quires p ^ and the dominant energy condition requires 
\Pr/p\ ^ lj \pe/p\ < 1 and \pz/p\ < 1. The parameters 
used in this article satisfy all energy conditions. Further- 
more, the level curves show that it is physically acceptable. 
We remark that these are not the only parameters in which 
the level curves are physically acceptable. In the next sec- 
tion we apply the selected perturbations of Sec. [2] to the 
general relativistic Miyamoto-Nagai disk mentioned above 
and study its stability. 



4 PERTURBATIONS 

Before applying the different kinds of perturbations to the 
general relativistic Miyamoto-Nagai disk we must do some 
considerations. Note that the general relativistic Miyamoto- 
Nagai disk is infinite in the radial and axial directions. We 
want to study the stability of a finite disk. So, in order to 
achieve this requirement we need a cutoff in the radial co- 
ordinate. In Eqs. (fl7)l . (jTSj and (|T9j , we see that the ther- 
modynamic quantities decrease rapidly enough to define a 
cutoff in both coordinates. The radial cutoff R cu t and the 
axial cutoff Z cu t are set by the following criterion: the en- 
ergy density within the disk formed by the cutoff parameters 
has to be more than 90% of the infinite thick disk energy 
density. The above criterion, and the parameters used in the 
article, leads to a radial cutoff of R cu t = 10 units and an ax- 
ial cutoff of \Z cut \ — 5 units. The other 10% of the energy 
density that is distributed from outside the cutoff parame- 
ters to infinity can be treated, if necessary, as a perturbation 
in the outermost boundary condition. 



4.1 Perturbation with SU, 5p e , SU R , 5p R 

We start perturbing the four velocity in its components 6 
and R. From the physical considerations mentioned in Sec. [2] 
we also expect variations in the thermodynamic quantities 
pe and pr. The set of equations I©-© reduces to a sec- 
ond order ordinary differential equation for the perturbation 
Sp R , say 



FaSpr.rr + F B 5p R R + F c SpR = 0, 



(20) 



where (Fa, Fb , Fc) are functions of (R, z,w, ke), see Ap- 
pendix fS] For this particular case we have, Spe = —SpR. 

Note that in Eq. (|20[) the coordinate z only enters as a 
parameter. Moreover, the equation for SpR is independent 
of the parameter w, but w needs to be different from zero 
to reach that form. The second order equation (|20|) is solved 
numerically with two boundary conditions, one at R ~ 
and the other at the radial cutoff. At R w we set the per- 
turbation SpR to be m 10% of the unperturbed pressure pr 
(|18|) . In the outer radius of the disk we set 5pR.\R=R cut = 
because we want our perturbation to vanish when approach- 
ing the edge of the disk and, in that way, to be in accordance 
with the applied linear perturbation. We say that our per- 
turbations are valid if their values are lower, or of the same 
order of magnitude, than the 10% values of its unperturbed 
quantities. 

In Fig. [T] we present the amplitude profile of the radial 
pressure perturbation in the plane z = for different values 
of the parameters a and b. As in the Newtonian case, the 
less the ratio b/a, the flatter is the mass distribution. We 
see that the perturbation SpR for (a — 1,6 = 1) decreases 
rapidly with R and has oscillatory behavior. At first sight, 
the perturbation SpR appears to be stable for all R, but in or- 
der to make a complete analysis we have to compare at each 
radius the values of the perturbations with the values of the 
radial pressure. For this purpose, we included in the same 
graph a profile of the 10% value of pr. We see that the per- 
turbations of SpR for different values of k g are always lower 
or, at least, of the same order of magnitude when compared 
to these 10% values. In the flatter case (a = 1,6 = 0.5), the 
perturbation SpR shows the same qualitative behavior, but 
the amplitudes of the oscillations are slightly higher. In both 
cases the amplitudes are well below the 10% values of pr. If 
we consider a very flat galaxy (a — 1,6 = 0.1) with w = 1 
we found that some modes are not stable in a small region 
near the center of the disk, from R ~ to R ~ 0.3, because 
the perturbation amplitude is bigger than the 10% value of 
Pr and our general linear perturbation is no longer valid. 

We also performed stability analyses for the physical ra- 
dial velocity perturbation SU = ^/gRRSU R and the phys- 
ical azimuthal velocity perturbation SU = ^/geeSU 6 . Note 
that our four velocity ?7 M @ has only components in the 
temporal part, so we do not have values of U R and U e to 
make comparisons with the perturbed values SU and SU . 
For that reason we compared, in first approximation, the 
amplitude profiles of these perturbations with the value of 
the escape velocity in the Newtonian limit. In the Newto- 
nian limit of General Relativity, / <C 1, we have the well 
known relation goo = ??oo + 2<E>. So, the Newtonian escape 
velocity V esc = W 2|$| can be written as V e sc = 2\/7, see 

IVogt &i Letelier! f2005bl ) . With this criterion, the perturba- 

~ r ~ e 

tions SU and SU are stable because their values are always 
well below the escape velocity value. Recall that the pertur- 
bation SpR does not depend on the parameter w, but the 
perturbations SU and SU do. We performed numerical 
solutions for the perturbations SU and SU with different 
values of the frequency w, and we find that when we increase 
the value of w the perturbations become more stable. 

In this subsection we set the value of the parameter 
z = 0. We performed the same analysis for different val- 
ues of the parameter —5 < z < 5, and we found that the 
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Figure 1. Profiles of the amplitude perturbation for the radial pressure of the fluid for the cases when z = 0, w = 1 and kg = 0, 1, 2, 3, 4. 
The graph at the left and right correspond to the cases when (a = 1, b = 1) and (a = 1, b = 0.5) respectively. The 10% p R profile is 
depicted for better stability comparisons. The modes of these examples are all stable. 



perturbations show the same qualitative behavior. There- 
fore, we can say that the general relativistic Miyamoto-Nagai 
disk shows some not-stable modes for very flat galaxies, e.g. 
(a — l,b = 0.1). Otherwise the disk is stable under pertur- 
bations of the form presented in this subsection. 

Nevertheless, if we treat the 10% of the energy density 
as a perturbation in the outermost radius of the disk by 
setting Sp R \ R=Rcut = e, where e < 10% of p R \ R=Rcut , the 
qualitative behavior of the mode profiles is the same. In the 
case of fiat galaxies, when they present not stable modes, 
more complex structures like rings, bars or halos can be 
formed. Moreover, if we set the frequency w — > iw we obtain 
the same equation for the perturbation Sp R , say (|20[) . In 
this case, the real part of the general perturbation diverges 
with time and the perturbation is not stable. These last 
considerations can be applied to every perturbation in the 
following subsections. 



4.2 Perturbation with SU , 5p e , SU Z , 5p z 

In this subsection we perturb the four velocity in its com- 
ponents and z, and we expect variations in the thermo- 
dynamic quantities pg and p z . The set of equations JS])-© 
reduces to a second order ordinary differential equation for 
the perturbation Sp z given by 

F A 8p z>zz + F B 8p Zz + F c 8p z = Q, (21) 

where (Fa, Fb , Fc) are functions of (R,z,w,kg), see Ap- 
pendix[B] Note that in Eq. (|2ip the coordinate R only enters 
as a parameter. Like the previous case, Eq. (|21fl is indepen- 
dent of the parameter w, but in order to reach that form we 
must have w different from zero. The second order equation 
()21[) is solved numerically with two boundary conditions, 
one in z — and the other in z — Z cu t. At z = we set the 
perturbation Sp z to be w 10% of the unperturbed pressure 
p z HHJ). In the outer plane of the disk we set Sp z \ z= z cut = 
because we want our perturbation to vanish when approach- 



ing the edge of the disk, and in that way, to be in accordance 
with the linear perturbation applied. 

In Fig [21 we present the amplitude profiles of the axial 
pressure perturbation, the physical axial velocity perturba- 
tion SU = y/g zz 8U z and the physical azimuthal velocity 
perturbation for R = 0.1 and different values of the pa- 
rameters a and b. For comparison reasons, we included in 
the graphs the amplitude profile that corresponds to 10% 
of the value of p z and the escape velocity profile. Note that 
for (a — 1, b — 1) some modes of the axial pressure per- 
turbation are above the 10% profile of p z , e.g. the modes 
with kg = and kg = 5. In these cases we can say that 
the mode with kg — is not stable and that the mode with 
kg — 5 is near the validity criterion used for the pertur- 
bations. These modes are also present in the flatter galaxy 
(a — 1,6 = 0.5) and have the same behaviors. The mode 
kg — 5 is actually not stable. This can be seen in the az- 
imuthal velocity perturbation profiles, where its amplitude 
is greater than the escape velocity. Note that in the velocity 
perturbation graphs the mode kg — is also not stable. The 
azimuthal pressure perturbation, not depicted in Fig [5] has 
all the modes well below the 10% profile of pg, and therefore 
is stable. 

The perturbations Sp z and Spg do not depend on the 
parameter w, but the perturbations SU and SU do. We 
performed numerical solutions for the perturbations SU 
and SU with different values of the frequency w, and we 
find that when we increase the value of w the perturbations 
become more stable. 

We have performed the same above analysis for different 
values of the parameter < R < 10, and we found that the 
qualitative behavior is the same. We see from Fig. [5] that 
the not stable modes are more pronounced for the flatter 
galaxy. Furthermore, for very flat galaxies some modes like 
kg = 10 become not stable. In general, for not stable modes, 
more complex structures like rings, bars or spiral arms may 
be formed. 
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Figure 2. Profiles of the amplitude perturbation for the axial pressure, axial physical velocity and azimuthal physical velocity of the 
fluid for the cases when R = 0.1, and kg = 0, 5, 10, 15, 20. The graphs at the left and right correspond to the cases when (a = 1, 6 = 1) 
and (a = 1,6 = 0.5) respectively. The 10% p z profile and the escape velocity are depicted for better stability comparisons. Note, in 
the axial pressure perturbation graph, that some nodes are not stable because they do not satisfy our stability criterion. In order to be 
stable, a mode must have the correct behavior in all the perturbed quantities. For example, in the case (a=l,b=l) the mode with w = 1 
and kg = 5 seems to be stable in 8p z , but looking into the SU 9 perturbation graph wc note that this statement is not true. In the natter 
galaxies the modes are not stable on larger regions of the domain. 
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Figure 3. Profiles of the amplitude perturbation for the radial pressure of the fluid for the cases when z = and w = 1, 5, 10, 15,20. 
The graphs at the left and right correspond to the cases when (a = 1,6 = 1) and (a = 1,6 = 0.5) respectively. The 10% Pr profile is 
depicted for better stability comparisons. These modes are highly stable. 



4.3 Perturbation with 5U R , 5pu, 5p 

In this subsection we perturb the radial component of the 
four velocity, the radial pressure and the energy density of 
the fluid. The set of equations ([5|-([8| reduces to a second 
order ordinary differential equation for the perturbation SpR 
of the form ((20j). The forms of the functions {F A ,F B ,F C ) 
are given in Appendix[C] In this case, the coordinate z only 
enters as a parameter. Due to the fact that we are not con- 
sidering perturbations in the azimuthal axis, the coefficients 
of the second order ordinary differential equation do not de- 
pend on the wavenumber kg. This second order equation is 
solved numerically with the same boundary conditions de- 
scribed in Sec. 14.11 

In Fig. Owe present the amplitude profiles for different 
perturbation modes of the radial pressure in the plane z = 
for different values of the parameters a and b. We see in the 
graph that the perturbation profiles decrease rapidly in few 
units of R. Also, the values of the radial velocity pertur- 
bation and energy density perturbation, not depicted, are 
well below the escape velocity and the 10% energy profile, 
respectively. 

We performed the above analysis for different values of 
— 5 < 2 < 5 and we found that the quantities involved have 
the same qualitative behavior. From these results, we can 
say that the general linear perturbation applied is highly 
stable and, for that reason, the perturbations do not form 
more complex structures. 



4.4 Perturbation with 5U Z , Sp z , Sp 

In this subsection we perturb the axial component of the four 
velocity, the axial component of the pressure and the energy 
density of the fluid. The set of equations (©-(HI reduces to 
a second order ordinary differential equation for the pertur- 
bation Sp z of the form (T2TJ). The functions (Fa,Fb,Fc) are 
given in Appendix [D] Note that, like in Sec. 14.21 the coor- 
dinate R only enters as a parameter. In this case, we are 



not considering azimuthal perturbations and therefore the 
quantities involved do not depend on the parameter kg. The 
second order equation is solved following the procedure of 
Sec. g3] 

In Fig. [4] we present the amplitude profiles of the axial 
pressure perturbation and the physical axial velocity pertur- 
bation, for R — 0.1 and for different values of the parameters 
a and b. We see that the axial pressure perturbation modes 
for (a = 1, b — 1) are always of the some order of magnitude 
or lower when compared to the 10% profile. In the flatter 
case (a = 1, b = 0.5), note that the amplitude of the mode 
10 = 1 is greater in some region of the domain. This fact 
is reflected in the axial velocity perturbation profile where 
the mode w = 1 have a strange behavior. All of the modes, 
including the mode with w = 1, are stable because they 
are well below the escape velocity, which is not depicted. 
The modes that correspond to the energy density perturba- 
tion are all stable. For highly flat galaxies the mode w — 1 
is not stable and may form more complex structures. For 
higher values of the parameter w the modes are more sta- 
ble. We performed the above analysis for different values of 
the parameter < R < 10 and we found that the quantities 
involved have the same qualitative behavior. 



4.5 Perturbation with 5U , SpR, SU Z , Sp z and 
8p R = Sp z 

In this subsection we perturb the radial component of the 
four velocity, the axial component of the four velocity, the 
radial pressure and the axial pressure. As we said in Sec. [2j 
we need an extra condition to set the number of unknowns 
equal to the number of equations. In this case, we set SpR = 
5p z = Sp. Therefore, the set of equations (©-(HI reduces to 
a second order partial differential equation for the pressure 
perturbation Sp, say 

F A Sp, RR + F B Sp, R + F c 5p, zz + F D Sp, z + F E 5p = 0, (22) 
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Figure 4. Profiles of the amplitude perturbation for the axial pressure and the axial physical velocity for the case when R = 0.1 and 
w = 1, 5, 10, 15, 20. The graphs at the left and right correspond to the cases when (a = 1, b = 1) and (a = 1, b = 0.5) respectively. The 
10% p z profile is depicted in the graph of 8p z for better stability comparisons. In the graph of 8U the escape velocity is not depicted 
because it is several orders of magnitude greater. For these examples all the modes are stable, but for a highly flat galaxy some modes, 
like the mode with u> = 1, are not stable. 



where (Fa, Fb , Fc , Fd , Fb) are functions of (R,z,w), see 
Appendix [E] The partial differential equation (|22|l is solved 
numerically with four boundary conditions, at z = —Z cu t, 
z — Zcut, R « and R — R C ut- They are different ways 
in which we can set the boundary conditions in order to 
simulate various kinds of pressure perturbations. Here, we 
treat only the case when we have a pressure perturbation at 
R « and along the z axis, i.e. some kind of a rod pertur- 
bation. We set the value of the rod pressure perturbation to 
be 10% of the axial pressure. We set the values of the other 
boundary conditions equal to zero because we want the per- 
turbation to vanish when approaching the edge of the disk. 
We choose the 10% of the value of the axial pressure instead 
of the radial pressure because it has the lowest value near 
R « 0. In that way, the perturbation values are also below 
the 10% values of the radial pressure and the general linear 
perturbation is valid. 



In Fig.[S] we present the perturbation amplitudes for the 
pressure, the physical radial velocity and the physical axial 
velocity, for w = 1 and for different values of the parame- 
ters a and b. We see in the pressure perturbation graph that 
the perturbation rapidly decays to values near zero when we 
move out from the center of the disk. This behavior is the 
same for every galaxy considered. In the velocity perturba- 
tions profiles we can see a phenomenon that is more clear in 
the flatter galaxy. Note that in the lower domain of the disk 
[-5,0) the axial velocity perturbation is positive and in the 
upper domain (0,5] the axial velocity perturbation is nega- 
tive. This means that due to the linear perturbation the disk 
tries to collapse to the plane z = 0. Now, if we look to the 
radial velocity perturbation graph, we note that the upper 
and lower parts depart from the center of the disk due to the 
positive radial perturbation. So, with these considerations, 
we may say that the disk tends to form some kind of ring 
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Pressure Perturbation Pressure Perturbation 




Figure 5. Profiles of the amplitude perturbation for the pressure, the radial physical velocity and the axial physical velocity of the fluid 
for the case when w = 1 and for a rod perturbation in R tss 0. With this kind of perturbation the disk tends to form some kind of ring 
around the center of the disk. This phenomenon is greater for highly flat galaxies and lower for more spherical systems. 
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around the center of the disk. This phenomenon is greater 
for highly fiat galaxies and lower for more spherical systems. 



5 CONCLUSIONS 

In this article we studied the stability of the re- 
cently proposed general relativistic Miyamoto-Nagai model 
[Vogt fc LeteTIei] (|2005bl )] by applying a general first order 
perturbation. We can say that the stability analysis per- 
formed is more complete than the stability analysis of par- 
ticle motion along geodesies because we have taken into ac- 
count the collective behavior of the particles. However, this 
analysis can be said to be incomplete because the energy 
momentum perturbation tensor of the fluid is treated as a 
test fluid and does not alter the background metric. This is 
a second degree of approximation to the stability problem in 
which the emission of gravitational radiation is considered. 

The different stability analyses made to the general rel- 
ativistic Miyamoto-Nagai disk show that this disk is stable 
for higher values of the wave number kg and the frequency 
w. For lower values of kg and w the disk presents not-stable 
modes that may form more complex structures like rings, 
bars or halos, but in order to study them we need a higher 
order perturbation formalism. In general, not-stable modes 
appear more for flatter galaxies and less for spherical sys- 
tems. 



APPENDIX B: FUNCTIONS F A , F B AND F c OF 
SECTION 4.2 

The general form of the functions (Fa, Fb,Fc) appearing in 
the second order ordinary differential equation (121 1) is given 
by 



FA = A 2 cti, F B = A 2 ai :Z + A 2 a 2 + Asai, 
F c = A 2 ct 2 ,z + A 4 a 3 + A 5 a 2 , 

where ai, a 2 and a 3 are 



(Bl) 



U 3 



_Di _ B 6 D S - B 5 D 6 

D 2 ' Q2 ~ B 5 D 2 
C 2 Bq 

els? 



(B2) 



and the meaning of the coefficients (Ai, Bi, C\, Di) is ex- 
plained in Appendix lAl 



APPENDIX C: FUNCTIONS F A , F B AND F c OF 
SECTION 4.3 

The general form of the functions (Fa, Fb, Fc) is given by 
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APPENDIX A: FUNCTIONS F A , F B AND F c OF 
SECTION 4.1 

The general form of the functions (Fa, Fb,Fc) appearing in 
the second order ordinary differential equation (f20|) is given 
by 



F A = A 1 a 1 , Fb = Aioti )R + Aia 2 + A 3 ai, 
F c = -4i«2,.r + A 3 a 2 + A 6 a 3 , 



(CI) 



where ai, a 2 and a 3 are 



ai = 



a 3 = 



Bi 
B 2 ' 
Fh 
D 3 



«2 = 



B 3 D 4 - B 4 D 3 
B 2 D 3 



(C2) 



and the meaning of the coefficients (Ai, Bi,Di) is explained 
in Appendix lAl 



F A =A 1 a 1 , F B = Aia 1M + A±a 2 + A 3 ai, 

F c = A!a 2 , R + A 3 a 2 + A 4 a 3 , (AI) 

where a±, a 2 and a 3 are 



APPENDIX D: FUNCTIONS F A , F B AND F c OF 
SECTION 4.4 

The general form of the functions (Fa, Fb, Fc) is given by 



Ql = - 



q 3 = 



B 2 ' 
C 2 Di 
C1D5' 



OL 2 



B 5 D 4 - BiD 5 
B 2 D 5 



(A2) 



In Eqs. (IA1|) and (|A2|) . we denote the coefficients of Eq. ([5]) 
by Ai, the coefficient of Eq. (|6]) by Bi, the coefficient of Eq. 
by d, the coefficient of Eq. ([HJ by A, e.g., the first 
term in <(Sj has the coefficient Ai multiplied by the factor 
SUji, the second term has the coefficient A 2 multiplied by 
the factor SU R , etc. The explicit form of the above equations 
is obtained replacing the fluid variables (p,pR,pg,p z ) of the 
isotropic Schwarzschild thick disk. 



F A = A 2 u\, F b = A 2 ai yZ + A 2 a 2 + A 5 ai, 
Fc = A 2 a 2 ^ z + Asa 2 + ^603, 

where oil, a 2 and a 3 are 



(Dl) 



B G D 3 - B 3 D 6 
a i — — -— , a 2 — — — . 

D 2 ' B 3 D 2 



a 3 



El 
D 2 

Be 

B 3 ' 



(D2) 



and the meaning of the coefficients (Ai, Bi,Di) is explained 
in Appendix \K\ 
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APPENDIX E: FUNCTIONS F A , F B , F c , F D 
AND F E OF SECTION 4.5 

The general form of the functions (Fa, Fb, Fc,Fd,Fe) ap- 
pearing in the partial second order differential equation (|22[) 
is given by 

Fa = Aiai, Fb = Aiai lR + Aiai + A3CX1, 
Fc — A2CX3, Fd = ^203,2 + ^204 + Ascw, 
F E = A ia2 , R + A 2 a 4)Z + A3CX2 + A 5 a t , (El) 

where ati, 012, 0.3 and Q4 are 



cti 



0t3 



Bi 

El 

D 2 > 



a 2 



Ct4 



B4 + Be 

B~2 ' 

D 2 ' 



(E2) 



and the meaning of the coefficients (Ai, Bi, Di) is explained 
in Appendix [A] 



REFERENCES 

Abramowicz MA., Prasanna A.R., MNRAS, 245, 720 
Bardeen J.M., 1970, ApJ, 161, 103 

Bicak J., Ledvinka T., 1993, Phys. Rev. Lett., 71, 1669 
Bicak J., Lynden-Bell D., Katz J., 1993, Phys. Rev. D, 47, 
4334 

Bicak J., Lynden-Bell D., Pichon C, 1993, MNRAS, 265, 
126 

Binney J., Tremaine S., 1987, Galactic Dynamics, Prince- 
ton University Press, New Jersey 

Bonnor W.B., Sackfield A., 1968, Commun. Math. Phys., 
8, 338 

Chazy J., 1924, Bull. Soc. Math. France, 52, 17 
Curzon H.E.J., 1924, Proc. London Math. Soc, 23, 477 
Frauendiener J., Klein C, 2001, Phys. Rev. D, 63, 084025 
Fridman A.M., Polyachenko V.L., 1984, Physics of Grav- 
itating Systems 1 : Equilibrium and Stability, Springer- 
Verlag, New York 
Garcfa-Reyes G., Gonzalez G.A., 2004, Phys. Rev. D, 69, 
124002 

Gonzalez G.A., Espitia OA., 2003, Phys. Rev. D, 68, 
104028 

Gonzalez G.A., Letelier P.S., 1999, Class. Quantum Grav., 
16, 479 

Gonzalez G.A., Letelier P.S., 2000, Phys. Rev. D, 62, 
064025 

Gonzalez G.A., Letelier P.S., 2004, Phys. Rev. D, 69, 
044013 

Kalnajs A.J., 1972, ApJ, 175, 63 

Karas V., Hure J-M., Semerak O., 2004, Class. Quantum 



Klein C, Richter O., 1999, Phys. Rev. Lett., 83, 2884 

Kuzmin G.G., 1956, Astrton. Zh., 33, 27 

Landau L.D., Lifshitz E.M., 1987, Fluid Mechanics, 2nd 

ed., Pergamon Press, Oxford, Sec. 27 
Ledvinka T., Zofka M., Bicak J., 1999, Proceedings of 

the 8th Marcel Grossman Meeting in General Relativity, 

edited by T. Piran, World Scientific, Singapore, 339 
Lemos J.P.S., 1989, Class. Quantum Grav., 6, 1219 
Lemos J.P.S., Letelier P.S., 1993, Class. Quantum Grav., 

10, L75 

Lemos J.P.S., Letelier P.S., 1994, Phys. Rev. D, 49, 5135 
Lemos J.P.S., Letelier PS., 1996, Int. J. Mod. Phys. D, 5, 

53 

Letelier PS., 1999, Phys. Rev. D, 60, 104042 

Letelier PS., 2003, Phys. Rev. D, 68, 104002 

Lord Rayleigh, 1916, Proc. R. Soc. London, A93, 148 

Lynden-Bell D., Pineault S., 1978, MNRAS, 185, 679 

Mestel L., 1963, MNRAS, 126, 553 

Miyamoto M., Nagai R., 1975, PAS J, 27, 533 

Morgan L., Morgan T., 1970, Phys Rev. D, 2, 2756 

Morgan T., Morgan L., 1969, Phys. Rev., 183, 1097 

Nagai R., Miyamoto M., 1976, PASJ, 28, 1 

Neugebauer G., Meinel R., 1995, Phys. Rev. Lett., 75, 3046 

Plummer H.C., 1911, MNRAS, 71, 460 

Satoh G., 1980, PASJ, 32, 41 

Seguin F.H., 1975, ApJ, 197, 745 

Semerak O., 2002, Gravitation: Following the Prague 
Inspiration, to Celebrate the 60th Birthday of Jiri 
Bicak, edited by Semerak O., Podolsky J., Zofka 
M., World Scientific, Singapure, 111, available at 
|http: / /xxx.lanl.gov/abs /gr-qc/0204025 1 
Semerak O., 2003, Class. Quantum Grav., 20, 1613 
Semerak O., Zacek M., 2000a, PASJ, 52, 1067 
Semerak O., Zacek M., 2000b, Class. Quantum Grav., 17, 
1613 

Toomre A., 1963, ApJ, 138, 385 

Ujevic M., Letelier P.S., 2004, Phys. Rev. D, 70, 084015 
Ujevic M., Letelier P.S., 2005, A&A, 442, 785 
Ujevic M., Letelier P.S., 2006, J. Comput. Phys., 215, 485 
Ujevic M., Letelier P.S., 2007, Gen. Relativ. Gravit., DOI: 

10.1007/sl0714-007-0438-y 
Vogt D., Letelier PS., 2003, Phys. Rev. D, 68, 084010 
Vogt D., Letelier P.S., 2004a, Class. Quantum Grav., 21, 

3369 

Vogt D., Letelier PS., 2004b, Phys. Rev. D, 70, 064003 
Vogt D., Letelier PS., 2005a, Phys. Rev. D, 71,084030 
Vogt D., Letelier PS., 2005b, MNRAS, 363, 268 
Vogt D., Letelier P.S., 2007, Triaxial Analytical Potential- 
Density Pair for Galaxies, PAJS, in press 
Voorhees B.H., 1970, Phys. Rev. D, 2, 2119 
Zipoy D.M., 1966, J. Math. Phys., 7, 1137 



Katz J., Bicak J., Lynden-Bell D., 1999, Class. Quantum 

Grav. 16, 4023 
King I.R., 1966, Astron. J., 71, 64 
Klein C, 1997, Class. Quantum Grav., 14, 2267 
Klein C, 2001, Phys. Rev. D, 63, 064033 
Klein C, 2002, Phys. Rev. D, 65, 084029 
Klein C, 2003a, Phys. Rev. D, 68, 027501 
Klein C, 2003b, Ann. Phys. (N.Y.), 12, 599 



